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Abstract 

The assumption of linear response of protein molecules 
to thermal noise or structural perturbations, such as 
ligand binding or detachment, is broadly used in the 
studies of protein dynamics. Conformational motions 
in proteins are traditionally analyzed in terms of nor- 
mal modes and experimental data on thermal fluctua- 
tions in such macromolecules is also usually interpreted 
in terms of the excitation of normal modes. We have 
chosen two important protein motors — myosin V and 
kinesin KIF1A — and performed numerical investiga- 
tions of their conformational relaxation properties within 
the coarse-grained elastic network approximation. We 
have found that the linearity assumption is deficient for 
ligand-induced conformational motions and can even be 
violated for characteristic thermal fluctuations. The de- 
ficiency is particularly pronounced in KIF1A where the 
normal mode description fails completely in describing 
functional mechanochemical motions. These results indi- 
cate that important assumptions of the theory of protein 
dynamics may need to be reconsidered. Neither a single 
normal mode, nor a superposition of such modes yield an 
approximation of strongly nonlinear dynamics. 



Author Summary 

Biological cells use a variety of molecular machines 
representing enzymes, ion channels or pumps, and mo- 
tors. Motor proteins are nanometer-size devices gener- 
ating forces and actively moving or rotating under the 
supply of chemical energy through ATP hydrolysis. They 
are crucial for many cell functions and promising for nan- 
otechnology of the future. Although such motors repre- 
sent single molecules, their operation cycles cannot be 
in detail followed in simulations even on the best mod- 
ern supercomputers and some approximations need to be 
employed. It is often assumed that conformational dy- 
namics of motor proteins is well described within a linear 
response approximation and corresponds to excitation of 
normal modes. We have checked this assumption for two 
motor proteins, myosin V and kinesin KIF1A. Our re- 
sults show that, while both these biomolecules respond 
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by well-defined motions to energetic excitations, these 
motions are essentially nonlinear. The effect is partic- 
ularly pronounced in KIF1A where relaxation proceeds 
through a sequence of qualitatively different conforma- 
tional changes, which may facilitate complex functional 
motions without additional control mechanisms. 



Introduction 

Protein machines, which may represent enzymes, ion 
pumps or molecular motors, play a fundamental role in 
biological cells and understanding of their activity is a 
major challenge. Operation of these machines is based on 
slow conformational motions powered by external energy 
supply, often with ligands (such as ATP). In molecular 
motors, binding of ATP and its subsequent hydrolysis in- 
duce functional mechanochemical motions, essential for 
their operation. These motions, which follow after an 
energetic activation, are conformational relaxation pro- 
cesses. 

Large-scale conformational changes may take place 
in proteins as a result of ligand binding [l]. Despite 
the large magnitude of such changes, they are nonethe- 
less often considered in the framework of the linear re- 
sponse theory [2[ and the normal mode approximation 
[3|-|7[. The normal mode analysis is furthermore broadly 
employed in the elastic-network studies of proteins 0- 
[Hj]. However, there is no general justification to as- 
sume that relaxation processes in proteins are linear and 
this assumption has to be verified for particular macro- 
molecules. 

It is known that relaxation processes in complex dy- 
namical systems may be strongly nonlinear and deviate 
much from simple exponential relaxation. As an exam- 
ple borrowed from a distant field, we can mention the 
Belousov-Zhabotinsky reaction which exhibits a great va- 
riety of spatiotemporal patterns (pacemakers, rotating 
spiral waves) that are however only complicated tran- 
sients accompanying relaxation to the equilibrium state 
[l9l [20| . Conformational relaxation in single protein 
molecules may also be a complicated process, compris- 
ing qualitatively different kinds of mechanochemical mo- 
tions. 

While partial unfolding and refolding, associated with 
ligand binding, are known for some protein machines, 
such as the enzyme adenylate kinase [21(, usually func- 
tional conformational motions in molecular machines 
and, specifically, in motor proteins are elastic. This 
means that the pattern of contacts between the residues 
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FIG. 1: Relaxation paths of myosin V. The elastic network model is constructed for the ATP-bound structure as the reference 
state. The red line shows the trajectory in the plane of distances 1112 and W13 between labels 1, 2 and 3 indicated in panel 
(A) starting from the nucleotide free state, so that this path corresponds to the conformational transition upon ATP binding. 
In panels (B) and (C) (magnified), gray lines display trajectories starting from 100 different random initial conditions (see 
Methods). In panel (D), gray lines represent relaxation trajectories with 100 different random deformations applied to the same 
initial structure as that for the red line. The dotted line in panel (D) shows the direction of distance changes corresponding to 
the slowest normal relaxation mode. Black dots indicate (meta)stable states reached. Times ti at points i — to 7 indicated 
in panel (C) are U = 0, 10, 30, 100, 300, 1000, 3000, and 10000, respectively. 




FIG. 2: Relaxation motion in myosin V. Snapshots of the conformations of myosin V along the special relaxation path are 
shown. The essential light chain, displayed in green, is included into the model. 
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in a protein is not changed upon ligand binding and 
preserved during the relaxation process, as generally as- 
sumed in the elastic network modeling (ENM) . 

Here, we provide detailed analysis of conformational 
relaxation processes, associated with ligand binding and 
hydrolysis, in two motor proteins — myosin V (22|, [23[ 
and kinesin KIF1A [24]. Our investigations, performed 
in the framework of the ENM approximation, reveal that 
nonlinearity is characteristic for both macromolecules 
and the normal mode description is not really applicable 
for any of them. For KIF1A, a monomeric motor protein 
from the kinesin superfamily, nonlinear effects are found 
to dominate completely functional mechanochemical mo- 
tions which turn out to be qualitatively different from the 
normal mode predictions. Despite the nonlinearity, well- 
defined conformational relaxation paths, robust against 
perturbations, have been found in both motor proteins. 



Results 

Within the coarse-grained ENM approach, a protein is 
modeled as a network of point-like particles, correspond- 
ing to residues, which are connected by a set of elas- 
tic links 0, [loj] ■ A link between two particles is present 
if the distance between them in the equilibrium confor- 
mation of the considered macromolecule is shorter than 
a cutoff length 

U = (n/2)Y Jij A ij (d ij -d, 
constant of the network links, Aaa is the matrix of con- 



The elastic energy of the network is 
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where k is the stiffness 



nections inside the network, is the distance between 

particles i and j, and is the respective distance in 
the equilibrium reference state. The characteristic time 
scales of functional mechanochemical motions in motor 
proteins are in the millisecond range and slow confor- 
mational relaxation motions on such timescale should 
be overdamped [25). Neglecting hydrodynamic interac- 
tions, relaxation dynamics is then described by equations 
dRi/dt = —TdU/dlli for the coordinates Ri of the par- 
ticles, where T is their mobility. Relaxation dynamics for 
clastic networks of proteins has been previously consid- 
ered 

Despite a wide-spread misunderstanding, elastic dy- 
namics is generally nonlinear. For example, macroscopic 
objects, such as ribbons or membranes, can still exhibit 
pronounced nonlinear effects of spontaneous twisting or 
buckling, while fully retaining their elastic behavior and 
not undergoing plastic deformations 27] . The energy U 
of an elastic network is quadratic in terms of the dis- 
tances dij and the forces acting on the particles are lin- 
ear in terms of such distances. However, the distance 
dij = |Rj — Rj| is itself a nonlinear function of the co- 
ordinates Ri and Rj and this makes the forces also non- 
linear functions of dynamical variables. The presence of 
nonlinear effects in conformational relaxation of proteins 
in the ENM approximation has been previously demon- 
strated [2ll2a. 



Explicitly, relaxation dynamics of considered proteins 
is described by equations (j3|) in the Methods section, 
where further details are also given. To study conforma- 
tional relaxation, these equations were numerically inte- 
grated starting from various initial conditions. 



Myosin V 

The reference conformation, used to construct the 
clastic network, was that of the ATP (analog) bound 
state (Protein Data Bank (PDB) ID code: 1W7J, with 
MgADP-BeFx as the ATP analog [30]). As the ini- 
tial condition, the conformation corresponding to the 
nucleotide-free state was taken (PDB ID: 10E9 |3lj). 
The elastic network had 855 particles connected by 7261 
links. Note that only the residues whose a-carbon po- 
sitions are contained in both PDB data sets have been 
taken to construct the network. Additionally, relaxation 
processes starting from randomly generated initial condi- 
tions (see Methods) have been considered. For visualiza- 
tion purposes, motions of three particles (Aspl22 in chain 
A, and Val22 and Serl35 in chain B) have been traced 
(Figure [JJA.) . Thus, each relaxation process was charac- 
terized by a certain trajectory in the space of distances 
between the three chosen labels. 

Figures [Tj3{Tp display 100 conformational relaxation 
trajectories, each starting from a different random initial 
condition. Although the initial conditions were gener- 
ated by applying relatively strong deformations (with- 
out unfolding) to the reference state, almost all of them 
were leading back to that reference state, with just a few 
metastable states found. Furthermore, one can observe 
that the trajectories converge to a well-defined relaxation 
path. 

The red trajectory in Figures [Tj3|Tp is for the relax- 
ation starting from the nucleotide-free conformational 
state of myosin V (so that the mechanochemical motion 
following ATP binding is simulated). After a transient, 
this special trajectory joins the well-defined relaxation 
path. This functional trajectory is robust against per- 
turbations, as shown by Figure [Tp. Several snapshots of 
the conformation along this trajectory are shown in Fig- 
ure^ (see also Video SI). 

The attractive path corresponds to a deep energy val- 
ley in the energy landscape of myosin V. Once this valley 
is entered, the conformational relaxation motion becomes 
effectively one-dimensional and characterized by a single 
mechanical coordinate. The profile of the elastic energy 
along the bottom of such energy valley determines the 
dependence of the elastic energy on the collective me- 
chanical coordinate (see Methods). 

Figure I3K shows the dependence of the elastic energy 
along the special attractive relaxation path starting from 
the nucleotide-free state and leading to the ATP-bound 
state. Markers indicate positions along the trajectory in 
Figure [Tp. For t > i 6 , the elastic energy U is approxi- 
mately quadratic in terms of the mechanical coordinate 
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FIG. 3: Profiles of elastic potential energy. (A) Myosin V and (B) KIF1A. Elastic energy U during transitions, corresponding 
to the trajectories shown by red lines in Figures [T] and [H is plotted against coordinate <f>. The dashed line shows U = ^Xi(j) 2 , 
the quadratic approximation corresponding to the slowest normal mode (see Methods). Numbers correspond to time moments 
indicated in Figures [Tp and[3p. 



</>, i.e. U(4>) ~ (A/2)0 2 . Because dcj)/dt = -dU/d<p, this 
implies that then d<fi/dt — — \<fi and the relaxation is ex- 
ponential. Only within such harmonical neighborhood of 
the reference state, the normal mode description becomes 
applicable (see Methods for further discussion). 

The dotted blue line in Figure [1)3 shows the direction 
of the distance changes corresponding to the slowest nor- 
mal mode (see Methods). The nucleotide-free state of 
myosin V lies away from this direction and also outside 
of the harmonical neighborhood. The initial stage of the 
functional mechanochemical motion (until time t§) can- 
not be quantitatively analyzed in terms of the normal 
modes. 



Kinesin KIF1A 

The reference conformation for KIF1A is the ADP- 
bound state (PDB ID: 1I5S, with MgADP [Hj]). Relax- 
ation starting from the initial condition, corresponding to 
the ATP(analog)-state (PDB ID: 1161, with MgAMPPCP 
as an ATP analog [32| ) and from randomly generated ini- 
tial conditions was considered. The elastic network has 
320 particles and 2871 links. Only the residues whose a- 
carbons are in both PDB data sets have been used. Visu- 
alization labels are Glu233, Ala286, and Asn211 (Figure 

m). 

100 relaxation trajectories, starting from random ini- 
tial conditions, are shown in Figures @}3 Up. The presence 
of an attractive relaxation path, corresponding to a deep 
energy valley, can be noticed. 

The red lines in Figures [4)3 |4p display the special re- 
laxation trajectory starting from the ATP-bound state. 
Surprisingly, we hnd that, in contrast to myosin V, this 
trajectory is different from the typical relaxation path. 



By applying small random initial perturbations to the 
initial ATP-bound state and integrating the dynamical 
equations, it can be demonstrated that this trajectory 
is, however, also stable with respect to the perturbations 
(FigureHp). The dotted blue line in Figure|3P shows the 
direction of the distance changes in the slowest normal 
mode of KIF1A. 

Thus, in KIF1A the deep energy valley leading to the 
reference ADP-bound state gets branched at some dis- 
tance from it. The path corresponding to the functional 
mechanochemical motion from the ATP-bound state be- 
longs to the side branch. Only at the final relaxation 
stage, in the immediate vicinity of the equilibrium, the 
valleys merge and the functional motion begins to coin- 
cide with the typical relaxation motion in this protein. 

The branching of the energy valley is already an indi- 
cation of strong nonlinearity in the relaxation dynamics. 
We have also determined the profile of the elastic energy 
U as a function of the mechanical coordinate <$> along the 
path connecting the ATP- and ADP-bound states (Fig- 
ure!^). The profile becomes quadratic only starting from 
time tg, very close to the equilibrium reference state. 

Figure [5] shows snapshots of KIF1A along the special 
attractive relaxation path (see also Video S2). At the 
early relaxation stage (until £5), the relaxation motion 
represents a combination of the rotation of the switch II 
helix and of the sliding of the switch I loop. Relaxation 
at the end of such initial stage is apparently hindered, as 
revealed in the presence of a plateau in the dependence of 
the elastic energy on the mechanical coordinate in Figure 
[3)3 near t = t$. Only once the sliding is completed, fur- 
ther local structural reorganization, representing a tran- 
sition from the loop to the cv-helix, becomes possible and 
is indeed observed approximately after time t§. 
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FIG. 4: Relaxation paths of KIF1A. The ADP-bound structure has been used to construct the elastic network. The visualization 
labels are indicated in panel (A), and the relaxation paths are displayed in panels (B) to (D) in the same way as in Figure Q] 
panels (B) to (D), respectively. The initial condition for the red line is the ATP-bound state, so that this trajectory corresponds 
to the transition from the ATP-bound state to the ADP-bound state. Times U at points i = to 9 indicated in panel (C) are 
U = 0, 0.1, 0.3, 1, 3, 10, 20, 25, 30, and 100, respectively. 



Discussion 



The normal mode description is broadly used in struc- 
tural studies of proteins. The analysis of thermal fluctua- 
tions and the interpretation of the respective experimen- 
tal structural data are traditionally performed assum- 
ing that fluctuations are linear and, hence, correspond 
to thermal excitation of various normal modes (see, e.g., 
[3, El)- The linear response of a protein macromolecule 
to structural perturbations, such as ligand binding, is an 
often used assumption [2|. To a large extent, the elastic- 
network analysis of ligand-induced macromolecular mo- 
tions is based on determining normal modes in the elastic 
networks of the considered proteins (see, e.g., [7]). The 
patterns of atomic displacements in such normal modes 
are further compared with the experimentally measured 
atomic displacements in the same proteins that are in- 
duced by a change of the chemical state, such as bind- 
ing of an ATP molecule 0, Il^ - [l7| . Large overlaps with 



only a few slowest normal modes are seen as the evi- 
dence for the applicability of the elastic-network ansatz, 
whereas the wide distribution of overlaps is considered as 
the indication that the elastic network description fails 
for a particular macromolecule. Specifically, strong over- 
laps between ligand-induced conformational changes and 
atomic displacements in the few slowest modes have been 
found for scallop myosin and Fi-ATPase, while such over- 
laps were absent for kinesin KIFIA [l8| . 

Our numerical investigations of elastic conformational 
motions in two motor proteins (myosin V and KIFIA) 
have revealed however that in both of them the nonlin- 
earities play an essential role. While slow conformational 
relaxation motions in myosin can still be qualitatively 
characterized in terms of the normal modes, the nor- 
mal mode description breaks down completely for KIFIA. 
The observed breakdown of the normal mode description 
does not however mean that conformational motions be- 
come irregular. We have seen that ordered and robust 
mechanochemical motions are characteristic for both pro- 
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FIG. 5: Relaxation motion in KIF1A. (A) Conformation snapshots (seen from two different viewpoints). Switch I and switch 
II regions are indicated in green and orange, respectively. (B,C) Schematic representation of the relaxation motion, observed 
in the simulation from the ATP-bound state (red) to the ADP-bound state (blue). In switch I region, as shown in panel (B), 
reconfiguration of the structure (2), i.e., transformation from a loop to an a-helix, occurs only after the sliding motion (1) is 
completed. For reference, relative positions of KIF1A and tubulin monomers (PDB ID: 2HXF and 2HXH [39|) are shown in 
panel (C). 
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FIG. 6: Eigenvalue spectra of the elastic network models. 
Eigenvalues X a in the normal mode description (eqn. (|10|l ). 
normalized to the lowest non-zero eigenvalue Ai, are shown. 
There is a significant gap between the lowest and the second 
lowest modes. 



tein motors, even though they cannot be described in 
terms of the linear response. 

We want to emphasize that, when the dynamics is non- 
linear, neither a single normal mode, nor a combination 
of many such modes can reproduce the motions. Thus, 
the normal mode description fails completely in this case 
and the problem is not that many normal modes must 
be taken into account. Actually, as we have shown, even 
for KIF1A, one normal mode would be sufficient to de- 
scribe long-time relaxation within the harmonic domain 
— however, this domain is restricted to a tiny neighbour- 
hood of the equilibrium state. 

Thermal fluctuations have not been explicitly included 
into our dynamical ENM simulations. However, such 
fluctuations are effectively generating random conforma- 
tional perturbations. In our study, relaxation processes 



starting from random conformational perturbations have 
indeed been considered. 

In myosin V, one well-defined nonlinear conformational 
relaxation trajectory, leading to the equilibrium state, 
has been identified. Starting from an arbitrary initial 
conformation (but still without unfolding), rapid conver- 
gence to this special trajectory takes place. While the 
motion corresponding to the special attractive trajectory 
is initially nonlinear, it becomes harmonical later and a 
substantial part of the ordered conformational relaxation 
process is within the harmonic domain of the equilibrium 
state. Similar behavior has been previously noted by us 
['29^ for scallop myosin and Fi-ATPase, but its detailed 
analysis has not yet been performed. 

The situation is more complex for the monomeric ki- 
nesin KIF1 A. Instead of a unique deep energy valley lead- 
ing to the reference ADP-bound state, two such valleys, 
both leading to the equilibrium state, arc present. These 
valleys correspond to two kinds of ordered conformational 
motions possible in the protein. 

The first of them is relatively wide and, when ther- 
mal conformational fluctuations are excited, they would 
typically proceed along it. However, the conformational 
relaxation motion starting from the ATP-bound state 
follows a different path, which corresponds to the sec- 
ond energy valley branching from the typical fluctua- 
tion path already at very small deviations from the equi- 
librium state. Note that the branching takes place as 
the change in the distance between the molecular labels 
Glu233 and Ala286 is still less than an angstrom, which 
is much smaller than the intensity of typical thermal flue- 
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tuations for such a distance. Thus, the nonlinear effects 
in KIF1A are strong even for the typical thermal fluctu- 
ations. 

Remarkably, such second relaxation path is also sta- 
ble with respect to perturbations, i.e. structurally ro- 
bust. Our numerical investigations reveal that motion 
along this path can be divided into two qualitatively dif- 
ferent stages. At the first of them, sliding of the switch 
I loop is observed, whereas at the second stage a tran- 
sition from the loop to the a-helix is realized. Struc- 
tural reorganization, corresponding to this transition, is 
not possible until the sliding motion is completed, lifting 
a restriction through the backbone chain. Recent crys- 
tallographic studies suggest that the switch I loop/helix 
plays an important role in control of the motor function 
through interaction with Mg 2+ and switch II [33[ • 

Thus, in contrast to myosin, a single ATP binding 
event induces in KIF1A a complex, but ordered confor- 
mational motion characterized by two qualitatively dif- 
ferent consequent phases. As we conjecture, this special 
dynamical property of KIF1 A may be needed for the pro - 
cessive motion of this single- headed molecular motor [2J] . 

In myosin V, conformational motions driven by ran- 
dom thermal fluctuations are similar in their proper- 
ties to the relaxation motion from the nucleotide-free 
state. This may facilitate exploitation of such fluctua- 
tion motions for the motor operation, as suggested by re- 
cent single- molecule experiments [34]. In KIF1A, where 
the energy valley splits into two branches, typical ther- 
mal conformational fluctuations are qualitatively differ- 
ent from the relaxation motion starting from the pre- 
vious ATP-bound state. The latter motion is entrop- 
ically hindered for thermal fluctuations and cannot be 
reversed through them. This may turn out to be im- 
portant for the understanding of the operation of the 
monomeric kinesin as a molecular motor. Latest exper- 
imental techniques permit simultaneous observation of 
stepping motion and conformational changes of a motor 
[351 ] . The coarse-grained modeling, including our present 
study, can contribute further suggestions for the design 
(e.g., by determining positions for fluorescent labeling) 
of such experiments. 

Finally, we note that our study has been based on 
the elastic network approximation for proteins. More 
detailed descriptions, such as, e.g., Go-like models, can 
also be used to consider conformational relaxation pro- 
cesses [36} . We expect that similar behavior will then be 
observed. 



the particles are determined by the locations of a-carbon 
atoms in the reference state of the protein, taken from 
the PDB database. Two particles in a network are con- 
nected by an elastic spring if at equilibrium the distance 



R 



(0) 



R 



(0) 



between them is less than a certain 



cutoff length Iq. The natural length of an elastic link is 
equal to the equilibrium distance d!f^ . The cutoff dis- 
tance Iq = 10 A has been used in our study. 

The elastic forces obey the Hooke law and all springs 
have the same stiffness constant n. Elastic torsion effects 
are not included. Thus, the force acting on particle i is 



N 



>(0) 



(i) 



where N is the total number of particles in the network, 
R-i(t) is the actual position of the particle i and dij = 
|R.j — Rj| is the actual distance between two particles i 
and j. The adjacency matrix of the network is defined as 
having — 1 if < Iq and Aij = otherwise. The 
total elastic energy of the network is 



(2) 



Because slow conformational dynamics of proteins in 
the solvent is considered, the motions are overdamped 
(see HU) and the velocity of a particle is proportional 
to the force acting on it, i.e. dRi/dt = TFi where T is 
the mobility. We assume that the mobilities of all parti- 
cles are the same. Hydrodynamical effects are neglected 
(they can be however incorporated into the elastic net- 
work models as shown in ref. [37|). 

Explicitly, the relaxation dynamics is described by a 
set of differential equations: 



N 



R, — R, 



R 



(0) 



R 



(0) 



R ? ; — Rj 

Rj — Rj | 

(3) 

Here, time is rescaled and measured in units of (Tk) . 
Hence, the relaxation dynamics of a network is com- 
pletely determined by its pattern of connections (matrix 
Aij) and the equilibrium distances d^ between the par- 
ticles. Equations (j3|) were integrated to determine con- 
formational relaxation motions. 



Methods 

Elastic network models 

In this study, we employ elastic network models where 
material points are connected by a set of elastic springs 
[§-[i~i"l]. Each particle corresponds to a residue in the 
considered protein. The equilibrium positions R^ of 



Simulations 

To prepare random initial conditions, the following 
procedure has been employed. Random static forces f^, 
acting on all particles in the network have been indepen- 
dently generated with the constraint that (1/N) ^ |fj | 2 = 
Ff ni - The equations of motion were integrated in the 
presence of such static forces for a fixed time r ini . The 
conformation which was thus reached has been then used 
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as the initial condition for the relaxation simulation. The 
parameters were F ini /K = 0.3 A, r ini — IC^ITk) -1 for 
myosin V and F ini /n = 0.5 A, T ini = 3 x 10 4 (Fk) _1 
for KIF1A. With these parameter values, relatively large 
overall deformations (~ 20% typical) could be reached, 
while still avoiding unfolding. In the deformed states, the 
lengths of the links did not exceed 13.4 A for myosin V 
and 11.5 A for KIF1A. 

When relaxation from specific conformations has been 
considered, initial positions of all particles were allocated 
according to the respective PDB structures. When ro- 
bustness of a relaxation path starting from a specific con- 
formation was investigated, the initial condition was pre- 
pared by randomly shifting the positions of all particles 
with respect to their locations in that conformation with 
a certain root- mean-square displacement We have 

chosen di n i = 2 A for myosin V and 0.5 A for KIF1A. 

To visualize conformational motions, three particles la- 
beled as 1, 2 and 3 were chosen and the distances u\i and 
Ui3 were monitored in the simulations. Thus, the relax- 
ation motion was represented by a trajectory on the plane 

(""12, W13). 

The choice of the visualization labels is essentially arbi- 
trary. In a simulation, motions of all residues were traced 
(see, e.g., Videos SI and S2) and different residues could 
be selected for a specific visualization. If a molecule has a 
low-dimensional attractive relaxation manifold, this is a 
property of the respective dynamical system and it can- 
not depend on the visualization method. When select- 
ing the labels, one should only pay attention to the fact 
that the distances between them should significantly vary 
during the relaxation process. If, by chance, two labels 
belonging to the same stiff domain in a protein have been 
taken, the distance between them would remain almost 
constant, so that such a choice would be inconvenient. 
When the normal mode description approximately holds 
and, furthermore, relaxation is well described by a few 
slowest modes, one can choose the labels so that the dis- 
tances between them reveal variations characteristic for 
the first few normal modes. Such selection was previ- 
ously made [29| for scallop myosin and Fi-ATPase, and 
it has been adopted in the present study for myosin V. 
For KIF1A, the labels have been chosen in such a way 
that motions in switch I and switch II regions are well 
resolved. 



Equation Q can be used to determine the coordinate cj> 
along a given relaxation trajectory and the dependence 
of the elastic energy U on this coordinate. 

For each point along the trajectory, the time t when 
it is reached in the relaxation process is known. More- 
over, the actual network configuration corresponding to 
this point is also known from the simulation. Therefore, 
for each point specified by time t the respective elastic 
energy U (t) is determined. The mechanical coordinate 
</>, reached at time t, is given by the integral 



0(t) 



dt 



(5) 



Normal mode description 

We provide a summary of the results on the normal 
mode description of conformational relaxation processes. 
If deviations = R, r' - 1 from the reference confor- 
mation are small for all particles, the nonlinear equations 
@, describing conformational relaxation of an elastic 
network, can be linearized: 



dt 



N 



,(°) 



u 



(0) 



(6) 



where u£ 0) = (bJ 0) - R?» ) /d^ 



Equations 



can be written in the matrix form as 



N 



^i--VAr- 



(7) 



where A is the 3N x 3A linearization matrix: 



A;, 



A i] U ij,a U ij 



v ( for i^j) 



V A .../CI 7, (0) 



(8) 



where a, 77 = (x, y, z). 

The general solution of these linear differential equa- 
tions is given by a superposition of 3N — 6 exponentially 
decaying normal modes, i.e. 



Profiles of elastic potential energy 



*(*) 



3N-6 
a=l 



(9) 



The collective mechanical coordinate <fi along a relax- 
ation path was defined by requiring that its dynamics 
obeys the equation d(f>/dt = —TdU/dcj) and that — > 
as t — > 00. Multiplying both parts of this equation by 
d(f>/dt, we find that it is equivalent to the equation 



dt 



dt ' 



(4) 



Here, A Q and are the eigenvalues and the eigenvec- 
tors of the linearization matrix, i.e. 



Ae 



(«) 



(a) 



(10) 



This matrix has 3N eigenvalues, but 6 of them must be 
zero, corresponding to free translations and rotations of 
the entire network. 
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Generally, all normal modes are initially present. As 
time goes on, first the normal modes with the larger 
eigenvalues A Q decay. In the long time limit, relaxation 
is characterized by the soft modes corresponding to the 
lowest eigenvalues. 

Figure [6] shows the computed eigenvalue spectra of 
myosin V and KIF1A. The eigenvalues are normalized 
to the lowest nonzero eigenvalue Ai and the logarithmic 
representation is chosen. 

Note that in both motor proteins a significant gap, sep- 
arating the soft mode from the rest of the spectrum, is 
present. This means that, in the linear approximation, 
long-time relaxation in these proteins is effectively char- 
acterized by a single degree of freedom, representing the 
amplitude of the soft mode. The pattern of displacements 
of particles (i.e., residues) from the reference positions is 
determined by the eigenvector of the soft mode. 

In the plane {u\%, U13) of the distances between the 
labels 1, 2 and 3, used by us for the visualization of 
conformational motions, the exponential relaxation mo- 
tion corresponding to the soft mode should proceed along 
the direction defined by the vector with the components 

u i2 ' (ei 1 ' ~ 4 X) ) and u i3 ■ { e i } ~ 4 1 ')- Such direc- 
tions are indicated by dotted blue lines in Figure QJ) and 
Figure Ht). 

When relaxation is reduced to a single soft mode, the 
elastic potential is quadratic in terms of the mechanical 
coordinate, i.e. U(4>) = (l/2)Ai</> 2 . 

Note that the representation of the relaxation process 
as a superposition @ of normal relaxation modes holds 
only in the harmonic domain, i.e. when linearization 
([6]) of full nonlinear relaxation dynamics equations ([3]) is 
valid. If dynamics is nonlinear and the linearization does 
not hold, relaxation dynamics cannot be viewed at all 
as a superposition of any normal modes. Whether just 
one normal mode or many of them should be included 
into a description of long-time relaxation dynamics is de- 
termined by the properties of the eigenvalue spectrum 
and not related to the possible invalidity of the harmonic 
approximation. 



As an extension,iterative normal mode analysis has 
been proposed [2l|, [38[ . This method is applied to obtain 
an optimal sequence of conformational states, transform- 
ing an initial given conformation into a target conforma- 
tion, which may be known with a low resolution or only 
partially, and thus to reconstruct missing details of that 
structure. Each next conformation in the sequence is ob- 
tained by making a step into the direction maximizing 
similarity with the target, restricted however to a super- 
position of a certain number of the lower normal modes. 
At the next iteration step, the previous conformation is 
chosen as a new reference state and a new set of normal 
modes is determined. This prediction method is useful 
and provides valuable results, e.g., in the refinement of 
low- resolution structures from electron microscopy [38l ]. 
It should be however emphasized that the sequence of 
conformational states yielded by such a method is gener- 
ally different from the path along which conformational 
relaxation from the target to the reference state would 
proceed. Even in the normal mode approximation, dy- 
namics of conformational relaxation depends not only on 
the eigenvectors, but also — and very significantly - 
on the eigenvalues of normal modes. Generally, the next 
iteration state in this method would not be the next con- 
formation along the actual relaxation path. This dif- 
ference can be clearly demonstrated by considering the 
example of KIF1A. The conformational relaxation path 
transforming the initial ATP-bound state into the (equi- 
librium) ADP-bound state is non- monotonous (Figure^) . 
It proceeds via intermediate states (particularly of the 
switch I region) which cannot be obtained by gradual in- 
terpolation maximizing similarity of the structures along 
the optimization path. 
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